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Recently the Yang-Mills gradient flow of pure SU(3) lattice gauge theory has been calculated in 
the range from P = 6/go = 6.3 to 7.5 (Asakawa et ah), where go is the bare coupling constant of the 
SU(3) Wilson action. Estimates of the deconfining phase transition are available from P = 5.7 to 6.8 
(Francis et al.). Here it is shown that the entire range from 5.7 to 7.5 is well described by a power 
series of the lattice spacing a times the lambda lattice mass scale An, using asymptotic scaling in the 
2-loop and 3-loop approximations for oKl- In both cases identical ratios for gradient flows versus 
deconfinement observables are obtained. Differences in the normalization constants with respect to 
A_l give a handle on their systematic errors. 

PACS numbers: 11.15.Ha 


I. INTRODUCTION 

We consider pure SU(N), N = 3, lattice gauge theory 
(LGT) with the Wilson action (see, e.g., Q) 

S = P = ^, ( 1 ) 

p ^0 

where the sum is over all plaquettes of a 4D hypercu- 
bic lattice with periodic boundary conditions. Up is the 
SU(N) plaquette variable, gg the bare coupling constant 
and P the usual convention, which emphasizes the in¬ 
terpretation as a 4D statistical mechanics, but gives up 
the P = l/{kT) relation with the physical temperature. 
Namely, T = l/{aNr) holds in LGT, where the integer 
Nr is the extension of the lattice in Euclidean time and 
a is the lattice spacing. 

For every physical observable m with the dimensions 
of a mass the relation 

m = Cm^■L (2) 

holds in the continuum limit a{P) -T 0 for /3 —>■ oo, where 
Ai sets the mass scale of the lattice regularization and 
Cm are calculable constants. Their actual computation 
faces difficulties, because one has to rely on simulations 
at finite lattice spacings a(/3), introducing corrections to 
the continuum relation. The subject of a good refer¬ 
ence scale arises. This topic gained renewed interest after 
Liischer Q introduced the Yang-Mills gradient flow scale, 
\/tQ, which comes by now in several variants. As antici¬ 
pated by Sommer in his review of the subject Q , gradient 
scales allow for an unprecedented precisiom when com¬ 
pared with traditional scales like ro or 0 defined by 
the force between static quarks at intermediate distance. 

In recent work Asakawa et al. @ pushed estimates for 
gradient scales in SU(3) gauge theory all the way up to 
P = 7.5. The SU(3) deconfining phase transition de¬ 
fines another precise scale, second only to gradient scales. 
Francis et al. @ managed to extend estimates of the 
SU(3) transition temperature Tt from lattice sizes of pre¬ 
viously Nr < 12 up to Nr = 22, Pt = 6.7986 (65). 


Remarkably, neither Asakawa et al. nor Francis et al. 
fit the P dependence of their estimates so that there is 
a /3 —> oo continuum limit as predicted by the universal 
part of asymptotic scaling. Instead, a parametrization 
for a limited P range is used and the continuum limit of 
ratios is subsequently estimated by fits in variables like 
(a/?'o)^) (a/\/to)^ and so on. This is in accord with a 
majority of publications on the subject, which all have 
given up on approaching the asymptotic scaling limit. 

Reasons for this, and why the decision to give up on 
asymptotic scaling may have been premature, are out¬ 
lined in section |IT1 Inspired by an earlier approach of 
Allton 0 , we are led to write the corrections to the mass 
relation ([2]) as a simple Taylor series in the lattice spacing 
times the lambda lattice mass scale, oA/,. In section uni 
this is seen to ;^ld excellent results for fitting the data 
of Ref. @ and p (see the abstract). Summary and con¬ 
clusions follow in the final section El 


II. ASYMPTOTIC SCALING AND 
CONTINUUM LIMIT 

The realization that the continuum limit of LGT 
may not just in theory but in practice be reached by 
computer simulations started with a paper by Creutz 
0, where he observed for the SU(2) string tension k 
a cross-over from its strong coupling behavior aPn = 
— ln(/3/4) to the I-loop asymptotic scaling behavior 
oPk = Ck exp(—67r^/3/ll). 

As the accuracy of Markov chain Monte Garlo calcu¬ 
lations improved, it was soon realized that there were, 
in particular for SU(3) with the Wilson action, strong 
violations of the asymptotic scaling relation and this did 
not improve noticeably by moving from the I-loop to the 
2-loop relation 

= fasisl) = exp : (3) 

where bo = llAI/(487r^) and bi = (34/3) AI^/(167r^)^ 
are, respectively, the universal I-loop I, IT^ and 2-loop 
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0[il coefficients of asymptotic freedom, called asymp¬ 
totic scaling in our context. Universal means that all 
renormalization schemes lead to the same bo and bi co¬ 
efficients. 

Next, the hope appeared to be that the situation would 
improve by including further, non-universal, terms of the 
expansion of aAj^: 


( 4 ) 

Computing up to 3-loops, Alles et al. [l^ calculated qi 
for SU(N) LGT, 

qi = 0.1896 forSU(3). (5) 

But, the discrepancies between the asymptotic scaling 
equation and data for physical quantities did not im¬ 
prove. 

Assuming that lattice artifacts are responsible for the 
disagreements, Allton Q suggested to include such cor¬ 
rections while constraining them with results from per¬ 
turbative expansions of the considered operators and ac¬ 
tions. Doubting, due to uncertainties with the very defi¬ 
nition of non-trivial continuum functional integrals, that 
perturbative information beyond Eq. o is reliable, a 
general Taylor series expansion in aA^, is proposed here 
for corrections to Eq. 


( OO 
1 + ^ 


m = Cm Al ^1 -I- 

where one has to determine the normalization constants 
Cm and the expansion coefficients di by computer sim¬ 
ulations. This has the potential to eliminate the essen¬ 
tial singularity of the perturbative expansion at = 0. 
However, the full sum (|3|) for /as(5o) available. In¬ 

stead, we have to work with approximations and define 
for g = 0, 1,... 


CO 


(oAl)* , oAl = /as(ffo) > (6) 


aAl = = f°s bo(/3)] 1 + ^ [9o(.P)y ) (7) 


i=i 


Expressing by /°j, gives rise to powers of logarithms 
like l/ln(/°^), ln|ln(/°j,)| and so on, which are singular 
for — >■ OO. Nevertheless, we continue to use with 

fas replaced by and come back to these issues after 
presenting the fits. 

In the following we consider observables with the di¬ 
mension of a length, L ^ 1/m and rewrite (jS]) as 


Lk 

a 


-1 


Cfe gAl I 1 + y^di (qAl)" 

V i=l 

C OO 

l + Ea.[/a.(5o")rj , 


( 9 ) 

( 10 ) 


where are the parameters with which we deal in our 
fits. There is no strong reason for using the expansion 
m instead of dH]). It just developed this way out of 
Ref. 0. To determine the expansion parameters Oi by 
numerical calculations one has to truncate the sum at 
rather small values of i. For sufficiently large /3 this 
should work well because {aAf) falls for all q exponen¬ 
tially off with /3 —>• OO. We define the truncated functions. 


iTW) = 




+Ear [fLm-^ , (11) 


with given by m and fit data according to 

= 'rw. ( 12 ) 

where the 2-loop {q = 0) and 3-loop (g = 1) asymptotic 
scaling functions, and ^j-e explicitly known 0. 
The labels p, q on the normalization constants Ck and 
parameters Ui indicate that their values depend on the 
choice of p, q. For simplicity the labels will be dropped 
when the association is obvious. 

For g = 0 as well as for g = 1 it turns out that excellent 
fits are obtained using p = 3 parameters at besides the 
Ck normalization constants. In the following we present 
g = 0,1, expansions for the Yang-Mills gradient flow 
data [1] and for the deconfining transition estimates Q. 


III. ANALYSIS OF THE NUMERICAL DATA 


where we have presently the g = 0 (2-loop) and g = I 
(3-loop) asymptotic scaling functions at our disposal 
and a conjecture for g 2 if we believe in the Bade approx¬ 
imation made in Ref. M- It is instructive to consider 
the deconfining temperature Tt as reference scale. Then 
a(A) = l/[N^{Pt)Tt] implies Ai(/3) = /«,(/3)iV,(/3)rt. 

Now, if the analyticity ® is true when using the full 
fas, it cannot be true at finite g. This is, for instance, 
seen by assuming that the expansion ([5]) is correct for 
f/g and comparing it with the same expansion using /°g. 
The difference lies in terms of the form 



(1 + 91 <70^-1 


( 8 ) 


For the gradient length scale a dimensionless variable 
f^{E{t)) is measured as a function of t. Then tx at 
which the observable takes a specific value X is used 
as reference scale. An operator whose t dependence 
has been extensively studied is E{t) = where 

Ffiv = — di,A^ -I- [A^, Ai/\ is the field strength. In 

Ref. @ solutions to the equations 


X and t^^t^{E{t)) 


A (13) 


have been calculated for X = 0.2, 0.3 and 0.4. The as¬ 
sociated l eng th scales are Vto. 2 i y/io.s, V'^qa and, intro¬ 
duced in [l5|, wo, 2 , wo. 3 , Wo A- For adaption to Eq. m 
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TABLE I: Error bars in percent of the signal, 100 ALfc/Lfc. 



Li 

L 2 

L 3 

L 4 

is 

Lq 


Lj 

p 

Vto.2 

Vto.3 

Vto.4 

Wo .2 

W 0.3 

Wo.4 

Pt 

Nr 

6.3 

0.09 

0.11 

0.12 

0.16 

0.17 

0.22 

5.69275 

0.07 

6.4 

0.07 

0.09 

0.08 

0.11 

0.12 

0.14 

5.89425 

0.05 

6.5 

0.13 

0.16 

0.19 

0.22 

0.21 

0.24 

6.06239 

0.06 

6.6 

0.12 

0.14 

0.16 

0.19 

0.21 

0.23 

6.20873 

0.07 

6.7 

0.26 

0.33 

0.35 

0.40 

0.46 

0.49 

6.33514 

0.06 

6.8 

0.18 

0.22 

0.25 

0.27 

0.30 

0.32 

6.4473 

0.25 

6.9 

0.46 

0.57 

0.65 

0.73 

0.81 

0.87 

6.5457 

0.54 

7.0 

0.14 

0.17 

0.19 

0.21 

0.25 

0.26 

6.6331 

0.26 

7.2 

0.43 

0.52 

0.59 

0.65 

0.71 

0.75 

6.7132 

0.34 

7.4 

0.30 

0.34 


0.41 

0.50 


6.7986 

0.84 

7.5 

0.37 



0.62 





rik 
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9 

11 

10 

9 
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they are renamed into Li,... ,Lq according to the first 
two rows of Table ID Their estimates are given in Table 1 
of @ and are not reproduced here. Instead, we give in our 
Table U error bars in percent of the signal, 100 
for the data tagged by a * in their paper, i.e., used in 
their analysis. 

Estimates of the SU(3) deconfining phase transition 
couplings Pt are given in Table I of Ref. @. Whenever 
(for smaller lattices) a comparison is possible their es¬ 
timates are consistent with previous work [H, [ 13 . The 
lengths associated with the deconfining phase transition 
temperatures Tt are l/(aTt) = Nt. However, the statis¬ 
tical errors are in Pt with Nr fixed. To allow for direct 
comparison with the other quantities, we attach to Nr 
error bars relying on the later estimated 1^^{P) scaling 
behavior from all data sets 


TABLE II: Xdof for our fits to each of the length scales. 


q 

Li 

L 2 

^3 

L 4 

As 

Lq 

Lr 

0 

0.46 

0.34 

0.23 

0.40 

0.47 

0.39 

0.76 

1 

0.42 

0.32 

0.24 

0.38 

0.46 

0.39 

0.74 


their error bars are systematically somewhat too large. 

For a visual presentation we have combined the entire 
n = ni -I- ... -|- ny = 70 data into two 1^'^{P), q = 0, 1, 
fits for Lk/{ack), which works astonishingly well. This 
is done with an extension of the method of [13. The 
constants Ck are defined as functions Cfe(ai, 02 , 03 ;data), 
which give the exact minimum of the fit for the particular 
constants Ui, effectively reducing the fitting procedure to 
three parameters, though the Ck are still counting against 
the degrees of freedom. The number of ai parameters is 
reduced by 6 x 3 to 3 from the 7 x 3 Oi parameters used 
altogether for the fits of Table [TTl 



P 


ANr = Nr 


lf(A + A/3,)-;f(A)l/;f(/30-(14) 


Starting with a guess and iterating the fit, one finds rapid 
convergence to the relative errors compiled in the Ly col¬ 
umn of Table m They are less than 0.25 for Pt < 6.33514 
{Nr < 12) and > 0.25 for Pt > 6.4473 {Nr = 14,..., 22), 
implying that the fit parameters will be dominated by 
the smaller Pt values. This is not good as the truncated 
parts of our expansion (HB) become more important at 
smaller p. Therefore, we adjust the Ly error bars for the 
lower Nr to 100 ALy/Ly = 0.2, which is still smaller than 
the best of the relative errors at the higher Nr values. 

For the gradient flow data the bias from smaller rela¬ 
tive errors is less severe and with /3 = 6.3 the smallest P 
is not so small. No adjustments are made in that case. 

The values of our fits (|T^ to the seven length 
scales are compiled in Table [IT] {ridof = n-k — 4: with Uk 
given in the last row of Table |I|. All fits are in very good 
agreement with the data. Actually, the fits of the gradi¬ 
ent flows are in too good agreement. This could be an 
accident, measurements of Li to Lq were performed on 
the same configurations so that they are all correlated, or 


FIG. 1: Data for I/y/(cya) and L 4 ,/{c 4 ,a) versus the fits 
(HU) using the 2-loop f^a (9 = 0) and the 3-loop /„„ {q — 1) 
asymptotic scaling functions. 


In Fig. [T] the two fits are shown jointly with the data 
points {i = 1,... ,nk) 

-^(*) ± ^^(f) for fc = 4, 7. (15) 

CL Cfi (X Cj^ 

Both fits cover with splendid X^of values the impressive 
range 5.69275 < P < 7.5. One value of k is picked for the 
gradient flow, because on the scale of the figure the data 
for the other Lk/{ack) lie right on top of them. For each 
q the first x^of value is for a fit that excludes the l/(aT 4 ) 
deconfinement data and the second x^^of value for the 
shown fit, which includes them. However, the increase 
from the values of Table [III should be noted. This 
and the fact that the data of Li to Lq are all correlated, 
as well as our “improvement” of the deconfinement data, 
may well obscure differences of the at parameters for dis¬ 
tinct observables. In fact, it is obvious from Fig. 4 (right) 



























4 


of Ref. @ that correlations greatly reduce the error bars 
of ratios and that \/t^/woA is not entirely flat as in our 
fits. To take these correlations into account one would 
best jackknife our fits, which requires the original time 
series. In the present context of simply demonstrating 
the almost identical scaling of all data graphically this 
would just be a distraction. Generally, one expects the 
oi parameters to agree for all so that corrections to 
ratios are of order (aA|^)^. There is no reason for 02 or 
03 to agree for all Lk- Only, it can be enforced within the 
accuracy of the present data. When these fits are applied 
to a single data set there is then a small bias due to the 
input of the other data sets. 

To make Fig. [T] reproducible, the fit parameters are 
given with high precision in Table IIIII More decent val¬ 
ues are obtained when one redefines the expansion pa¬ 
rameters aA|^ by multiplicative constants, e.g., so that 
they become 1 at /3 = 6 , = HsiP )/The 

second row of Table Ell gives the fits parameters for 
this case with their error bars in the third row. The 
range covered by the a;'^(/3) goes from x'^{5.7) « 1.4 
down to x'^{7.5) ~ 0.18, so that x'^{7.5)‘^ « 0.032 and 
a;'*(7.5)^ « 0.0058 become really small. 


TABLE III: Fit parameters used for Fig. [T] 








-155.559 

24615.3 

-5834850 

-104.735 

9926.28 

-2673493 

-0.365 

0.135 

-0.754 

-0.292 

0.773 

-0.581 

(13) 

(20) 

(82) 

(13) 

(42) 

(82) 


0.02 

0.015 

0.01 

0.005 

-A 0 

CO^ 

^ - 0.005 

<1 

- 0.01 

- 0.015 

- 0.02 

- 0.025 


FIG. 2: Relative deviations m of the L 7 and L 4 data points 
from the ht. 



5.8 6 6.2 6.4 6.6 6.8 7 7.2 7.4 




Fig. [ 2 ] provides a visual impression for the quality of 
the fits, by plotting the deviations of the fc = 4 and 7 
data points from the g = 1 fit of Fig. [T]in the form 


with A/cZ?’®(z) 


^(*)-zr(A) ( 16 ) 

ack 


together with error bars ALk/(ackl^'^)- 


TABLE IV: Normalization constants x 100. 


q 

0 

1 

0 

1 

Cl 

0.4918 (37) 

0.5569 (41) 

0.492 (06) 

0.557 (07) 

C2 

0.6198 (46) 

0.7018 (52) 

0.619 (11) 

0.701 (12) 

C 3 

0.7152 (53) 

0.8099 (59) 

0.695 (25) 

0.789 (28) 

C4 

0.5392 (40) 

0.6106 (45) 

0.548 (10) 

0.620 (11) 

C6 

0.6304 (47) 

0.7139 (52) 

0.638 (16) 

0.722 (17) 

C6 

0.7028 (53) 

0.7958 (59) 

0.679 (34) 

0.771 (38) 

C7 

2.4404 (71) 

2.7754 (79) 

2.357 (46) 

2.693 (51) 


Perhaps surprisingly, instead of one satisfactory de¬ 
scription of the data we got two (seven more pairs for 
the fits with their Xdof values listed in Table [TTI- The 
quality of the fits does not care about the log corrections 
discussed after Eq. ([ 8 ]). Instead, the parameters adjust 
and the normalization constants Ci to C 7 get shifted as 
shown in Table lYl Here the numbers in column 2 and 3 
correspond to the joint fits of the six gradient flow opera¬ 
tors, with exception of the last row, which corresponds to 
the fits displayed in Fig. [T] for which all seven operators 
are combined. Columns 4 and 5 give the results obtained 
from individual fits to which one should fall back when 
it comes to conservative estimates. Normalization con¬ 
stants of corresponding q = 0,1 fits differ by about 12 %, 
while their statistical errors are much smaller. 


TABLE V: Ratios of normalization constants. 


q 

0 

1 

0 

1 

C 7 

0.19728 (22) 

0.19724 (22) 

0.209 (05) 

0.207 (05) 

c? 

0.24861 (28) 

0.24856 (28) 

0.263 (07) 

0.260 (07) 

C 7 

0.28689 (33) 

0.28683 (33) 

0.295 (12) 

0.293 (12) 

C 7 

0.21630 (26) 

0.21625 (26) 

0.233 (07) 

0.230 (06) 

4 

0.25288 (32) 

0.25283 (32) 

0.271 (09) 

0.268 (09) 

4 

0.28188 (37) 

0.28182 (37) 

0.288 (16) 

0.286 (15) 


For ratios, cf = Ck/ci, of the normalization constants 
these differences become tiny and are swallowed by the 
statistical error bars as is seen in Table|V]for Cy (columns 
are arranged as in Table ITvT) . The deconfining transition 
is used as reference scale, because L 7 is statistically in¬ 
dependent from Li to Lq. The estimates of the last row 
can be compared with Asakawa et al. [^. Using q = 1, 
our values c® = woATt = 0.28182 (37) and 0.286 (15) are 
both well consistent with 0.285 (5) as given in their Ta¬ 
ble 3. Our value from column 3 is inconsistent with the 
precise estimate given in their Eq. (3.2), 0.2826 (3). The 
discrepancy may be well explained by the small bias of 
our result and/or the fact that Asakawa et al. rely en¬ 
tirely on Nr = 12 , whereas here a continuum fit is used 
that gives weight to all lattices, including Nr = 14 to 22. 

The Xdof of (fT^ are not sensitive to including 

or not including the qig^ term into the scaling function 
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FIG. 3: Approach of the fits to the asymptotic 2-loop and 
3-loop scaling functions 1°’'^ (q — 0,1) times 100/1°’^. 


while there is a remarkable shift in the normalization 
constants. It is then tempting, but entirely wrong, to 
argue that the gQ dependence is so weak that it does not 
matter and one could replace qigQ by a constant, say 
Qi9o Qi = 0-9 Qi for our /3 range. It is easy to see that 
with this, or any other q^, the normalization constants 
of the q = 0 fits will not change at all. So, the shift in 
the normalization constants comes entirely from the g^ 
dependence of the qi term. These contributions re-sum 
in a way that they become for large /3 responsible for the 
difference between and 

We use our fits of all n = 70 data to illuminate the situ¬ 
ation by Fig. [21 where for 5 = 0,1 the inverse asymptotic 
scaling functions and their fits are plotted times 
100/ i.e., as fractions of the inverse 3-loop asymptotic 

scaling function . We see that the gap between the 
and asymptotic scaling functions narrows slowly 
and the fits and approach rapidly (exponentially 
fast for increasing /3) their respective asymptotic behav¬ 
iors, where the fit stays closer to its asymptotic form 
than the fit: r; 0.8Z°’°/Z^’° over the entire /3 

range of the figure.. 

How does it come that the data cannot figure out 
whether the q = 0 or q = 1 fit is better? The answer lies 
in their ratios: If the ratio of the two fits is a constant, 
the difference between them will be entirely absorbed by 
the normalization. Defining the change in the ratios with 
respect to /3 = 6 as reference point by 


dP{l3) = 100 




(17) 


we find for the asymptotic scaling {p = 0 ) functions a 
change by 3.2% at /3 = 7.5. With 0.16% it is twenty 
times smaller for the fits (p = 3). 

What is then the effect of including more and more 
qj terms in the expansion (|T]) of /qs? We may expect 


convergence of the resulting normalization constants Ck 
towards their correct value. But how fast? Repeating 
the fits of all data with fake functions © defined by 
52 = ±0.19, so that 52 has a similar absolute value as 51 , 
there is again no sensitivity of the Xdof of the fits for the 
additional term and corrections to the Ck normalization 
constants stay less than ±10%. On this basis we end up 
with the result that our most reliable estimates of the 
Ck are those of column five of Table m with a mainly 
systematic uncertainty of ±10%. From c^ we get 

Al/Tt = C 7 ± 10% = 0.0269 (27) (18) 

in good agreement with Francis et al. who give 
Tt/Aj^g = 1.24(10). Using standard relations between 
lambda scales [T| this becomes A^/Tt = 0.0280 (25). Sim¬ 
ilarly, our estimate for wqaAl, 

LqA\ = C 6 ± 10% = 0.0077 (9), (19) 

is in agreement with the one of Table 3 of Asakawa et al. 
0 and the more accurate value of their Eq. (3.3), which 
translate, respectively, into wn, 4 AL = 0.00809(35) and 
woaAl = 0.00829 (5). 

When we believe in the Fade approximation of 
we find 52 = —0.02467, which is in magnitude almost ten 
times smaller than the range we allowed for our estimate 
of the systematic error. Using then fits with fasiP) as 
reference, Eq. m and m improve to 

Al/Tt = 0.0266 (9) and LqAI = 0.00762 (45), (20) 

where contributions of the statistical errors exceed now 
the systematic errors. So, it is difficult to understand why 
the error in Eq. 3.3 of Asakawa et al. is much smaller. 
Anyway, a small 52 suggests rapid convergence of the 
systematic errors of the normalization constants under 
increasing q for the used f^g functions. 


IV. SUMMARY AND CONCLUSIONS 

It appears that Eq. ([6]) is a natural parametrization 
of lattice spacing corrections to the continuum limit of 
SU(3) LGT. Incorporation of asymptotic scaling is still 
a viable alternative to other fitting methods for the ap¬ 
proach to the continuum limit, which are utilized in 
and elsewhere. In a next step, our fitting procedure 
should be tested for other asymptotically free theories, 
in particular full QCD. 


Acknowledgments 


This work was in part supported by the US Depart¬ 
ment of Energy under contract DE-FG02-13ER41942. I 
would like to thank David Glarke for calculating 52 from 
the Fade approximation of 0- 











6 


[1] I. Montvay and G. Munster, Quantum Fields on a Lat¬ 
tice, Cambridge University Press, 1994. 

[2] M. Liischer, JEEP 08, 071 (2010); 03, 092(E) 2014. 

[3] R. Sommer, POS (Lattice 2013) 015. 

[4] S. Necco and R. Sommer, Nucl. Phys. B 622 (2002) 328. 

[5] M. Asakawa, T. Hatsuda, T. Iritani, E. Itou, M. Ki- 
tazawa, and H. Suzuki, arXiv:1503.06516Y2. 

[6] A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus, and 
H. Ohno, Phys. Rev. D 91, 096002 (2015). 

[7] C.R. Allton, Nucl. Phys. B (Proc. Suppl.) 53, 867 (1997). 

[8] M. Creutz, Phys. Rev. D 21, 2308 (1980). 

[9] D.J. Gross and F. Wilczek, Phys. Rev. Lett. 30, 1343 
(1973). 

[10] H.D. Politzer, Phys. Rev. Lett. 30, 1346 (1973). 

[11] D.R.T. Jones, Nucl. Phys. B 75, 531 (1974). 


[12] W. Gaswell, Phys. Rev. Lett. 33, 244 (1974). 

[13] B. Alles, A. Feo and H. Panagopoulos, Nucl. Phys. B 491, 
498 (1997). 

[14] M. Gockeler, R. Horsley, A.C. Irving, D. Pleiter, P.E.L. 
Rakow, G. Schierholz, and H. Stiiben, Phys. Rev. D 73, 
014513 (2006). 

[15] S. Borsanyi, S. Diirr, Z. Fodor, C. Hoebling, S.D. Katz, 
S. Krieg, T. Kurth, L. Lellouch, T. Lippert, G. McNeile, 
and K.K. Szabo, JEEP 09, 010 (2012). 

[16] G. Boyd, J. Engels, F. Karsch, E. Laermann, G. Leg- 
eland, M. Liitgemeyer and B. Peterson, Nucl. Phys. B 
469, 419 (1996). 

[17] B.A. Berg and H. Wu, Phys. Rev. D 88, 074507 (2013). 

[18] B.A. Berg, arXiv: 1505.07564, submitted to CPG. 



